Take-Home Excercise 01

Author

Mayuri Salunke

Published

January 30, 2023

Modified

February 11, 2023

1.0 Overview

Water is an essential part of our life, without which we won’t be able to survive for more than 3 days. Living in Singapore, we have access to clean drinkable water 24/7 and as a result we don’t realize the struggles of people who don’t have access to clean water at all. An an example of such are the people in Nigeria. Despite 70% of Nigerians having access to basic water services, more than half of them are contamintated (Reference).

For this assignment, we will be focusing on the State of Osun in Nigeria. Osun, located in the southwestern Nigeria is bounded to the east by Ekiti and Ondo states, Kwara on the north, Ogun to the south and to the west by Oyo State (Reference). Their economy is mainly based on the agriculture and it inhibits the Osun River, a sacred river. However, in the recent years, the river has been polluted by the several mining activities from the surrounding communities (Reference). Hence, it is integral for us to address the issue of providing clean and sustainable water to the people of Osun. Through this assignment, I aim to apply the relevant spatial point pattern analysis learned in class to analyse the Functional and Non-Functional water points in State of Osun, Nigera.

Osun River

2.0 Setup

2.1 Packages Used

  • sf : Used for importing geospatial data, assigning or transforming coordinate systems, and converting geospatial and aspatial data into a sf data frame

  • tidyverse : Used for transforming and better presentation of Data

  • tmap : Used for plotting static point patterns maps or interactive maps

  • spatstat : Used for point-pattern analysis

  • raster : Used to read, write, manipulate, analyse and model gridded spatial data

  • maptools : Used to provide a set of tools for manipulating geographic data

  • rgdal : ??? DID I USE IT?

  • kableExtra : Used for table customization

  • funModeling : Used to data cleaning, importance variable, analysis and model performance

  • sfdep : FINDDD OUTTTT!!

pacman::p_load(sf, maptools, raster, spatstat, tmap, kableExtra, tidyverse, funModeling, sfdep)

2.2 Datasets Used

The below diagram shows the datasets used for the Assignment. We have two types of data - geospatial and aspatial.

For the Aspatial data, we are extracting the data from WPdx Global Data Repositories. The data source consists of two types of data - WPdx-Basic and WPdx+, for the purpose of this project, we will be using the WPdx+.

For the Geospatial data, we will be using the Nigeria Level-2 Administrative Boundary polygon features GIS data. There are two data source for this - Humanitarian Data Exchange (HDE) and geoBoundaries.

# initialise a dataframe of our geospatial and aspatial dataset details
datasets <- data.frame(
  Type=c("Geospatial",
         "Geospatial",
         "Geospatial",
         "Geospatial",
         "Geospatial",
         "Geospatial",
         "Geospatial",
         "Geospatial",
         "Geospatial",
         "Geospatial",
         "Geospatial",
         "Geospatial",
         "Geospatial",
         "Geospatial",
         
         "Aspatial"),
  
  Name=c("geoBoundaries-NGA-ADM2",
         "geoBoundaries-NGA-ADM2",
         "geoBoundaries-NGA-ADM2",
         "geoBoundaries-NGA-ADM2",
         "geoBoundaries-NGA-ADM2",
         "geoBoundaries-NGA-ADM2",
         "nga_admbnda_adm2_osgof_20190417",
         "nga_admbnda_adm2_osgof_20190417",
         "nga_admbnda_adm2_osgof_20190417",
         "nga_admbnda_adm2_osgof_20190417",
         "nga_admbnda_adm2_osgof_20190417",
         "nga_admbnda_adm2_osgof_20190417",
         "nga_admbnda_adm2_osgof_20190417",
         "nga_admbnda_adm2_osgof_20190417",
         
         "WPdx"),
  
  Format=c(".dbf", 
           ".geojson", 
           ".prj", 
           ".shp", 
           ".shx", 
           ".topojson",
           ".CPG",
           ".dbf",
           ".prj",
           ".sbn", 
           ".sbx", 
           ".shp", 
           ".shp", 
           ".shx", 
          
           ".csv"),
  
  Source=c("[geoBoundaries](https://www.geoboundaries.org/index.html#getdata)",
           "[geoBoundaries](https://www.geoboundaries.org/index.html#getdata)",
           "[geoBoundaries](https://www.geoboundaries.org/index.html#getdata)",
           "[geoBoundaries](https://www.geoboundaries.org/index.html#getdata)",
           "[geoBoundaries](https://www.geoboundaries.org/index.html#getdata)",
           "[geoBoundaries](https://www.geoboundaries.org/index.html#getdata)",
           
          "[Humanitarian Data Exchange](https://data.humdata.org/dataset/cod-ab-nga)",
           "[Humanitarian Data Exchange](https://data.humdata.org/dataset/cod-ab-nga)",
           "[Humanitarian Data Exchange](https://data.humdata.org/dataset/cod-ab-nga)",
           "[Humanitarian Data Exchange](https://data.humdata.org/dataset/cod-ab-nga)",
           "[Humanitarian Data Exchange](https://data.humdata.org/dataset/cod-ab-nga)",
           "[Humanitarian Data Exchange](https://data.humdata.org/dataset/cod-ab-nga)",
           "[Humanitarian Data Exchange](https://data.humdata.org/dataset/cod-ab-nga)",
           "[Humanitarian Data Exchange](https://data.humdata.org/dataset/cod-ab-nga)",
           
           "[ WPdx Global Data Repositories](https://www.waterpointdata.org/access-data/)")
  )

# with reference to this guide on kableExtra:
# https://cran.r-project.org/web/packages/kableExtra/vignettes/awesome_table_in_html.html
# kable_material is the name of the kable theme
# 'hover' for to highlight row when hovering, 'scale_down' to adjust table to fit page width
library(knitr)
library(kableExtra)
kable(datasets, caption="Datasets Used") %>%
  kable_material("hover", latex_options="scale_down")
Datasets Used
Type Name Format Source
Geospatial geoBoundaries-NGA-ADM2 .dbf [geoBoundaries](https://www.geoboundaries.org/index.html#getdata)
Geospatial geoBoundaries-NGA-ADM2 .geojson [geoBoundaries](https://www.geoboundaries.org/index.html#getdata)
Geospatial geoBoundaries-NGA-ADM2 .prj [geoBoundaries](https://www.geoboundaries.org/index.html#getdata)
Geospatial geoBoundaries-NGA-ADM2 .shp [geoBoundaries](https://www.geoboundaries.org/index.html#getdata)
Geospatial geoBoundaries-NGA-ADM2 .shx [geoBoundaries](https://www.geoboundaries.org/index.html#getdata)
Geospatial geoBoundaries-NGA-ADM2 .topojson [geoBoundaries](https://www.geoboundaries.org/index.html#getdata)
Geospatial nga_admbnda_adm2_osgof_20190417 .CPG [Humanitarian Data Exchange](https://data.humdata.org/dataset/cod-ab-nga)
Geospatial nga_admbnda_adm2_osgof_20190417 .dbf [Humanitarian Data Exchange](https://data.humdata.org/dataset/cod-ab-nga)
Geospatial nga_admbnda_adm2_osgof_20190417 .prj [Humanitarian Data Exchange](https://data.humdata.org/dataset/cod-ab-nga)
Geospatial nga_admbnda_adm2_osgof_20190417 .sbn [Humanitarian Data Exchange](https://data.humdata.org/dataset/cod-ab-nga)
Geospatial nga_admbnda_adm2_osgof_20190417 .sbx [Humanitarian Data Exchange](https://data.humdata.org/dataset/cod-ab-nga)
Geospatial nga_admbnda_adm2_osgof_20190417 .shp [Humanitarian Data Exchange](https://data.humdata.org/dataset/cod-ab-nga)
Geospatial nga_admbnda_adm2_osgof_20190417 .shp [Humanitarian Data Exchange](https://data.humdata.org/dataset/cod-ab-nga)
Geospatial nga_admbnda_adm2_osgof_20190417 .shx [Humanitarian Data Exchange](https://data.humdata.org/dataset/cod-ab-nga)
Aspatial WPdx .csv [ WPdx Global Data Repositories](https://www.waterpointdata.org/access-data/)

3.0 Data Wrangling : Geospatial Data

3.1 Importing and Transforming Geospatial Data

We will begin by importing Geospatial data into R by using the st_read() of sf package. It imports the nga_admbnda_adm2_osgof_20190417 shapefile into R as a polygon data frame. We provide 2 arguments - dsn (which is the data path) and layer (the shapefile name)

We use the st_transform() to perform projection transaction.

geoNGA <- st_read("data/geospatial", 
                  layer = "geoBoundaries-NGA-ADM2")
Reading layer `geoBoundaries-NGA-ADM2' from data source 
  `C:\mayurims\IS415-GAA\Take-Home_Ex\Take-Home_Ex01\data\geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 774 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 2.668534 ymin: 4.273007 xmax: 14.67882 ymax: 13.89442
Geodetic CRS:  WGS 84
NGA <- st_read(dsn = "data/geospatial", 
                  layer = "nga_admbnda_adm2_osgof_20190417")
Reading layer `nga_admbnda_adm2_osgof_20190417' from data source 
  `C:\mayurims\IS415-GAA\Take-Home_Ex\Take-Home_Ex01\data\geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 774 features and 16 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 2.668534 ymin: 4.273007 xmax: 14.67882 ymax: 13.89442
Geodetic CRS:  WGS 84

We can use the glimpse() of dplyr to know more about the associated attribute information of the dataframe.

glimpse(geoNGA)
Rows: 774
Columns: 7
$ shapeName  <chr> "Aba North", "Aba South", "Arochukwu", "Bende", "Ikwuano", …
$ pcode      <chr> "NG001001", "NG001002", "NG001003", "NG001004", "NG001005",…
$ level      <chr> "ADM2", "ADM2", "ADM2", "ADM2", "ADM2", "ADM2", "ADM2", "AD…
$ shapeID    <chr> "NGA-ADM2-13203401B25860527", "NGA-ADM2-13203401B76240303",…
$ shapeGroup <chr> "NGA", "NGA", "NGA", "NGA", "NGA", "NGA", "NGA", "NGA", "NG…
$ shapeType  <chr> "ADM2", "ADM2", "ADM2", "ADM2", "ADM2", "ADM2", "ADM2", "AD…
$ geometry   <MULTIPOLYGON [°]> MULTIPOLYGON (((7.387495 5...., MULTIPOLYGON (…
glimpse(NGA)
Rows: 774
Columns: 17
$ Shape_Leng <dbl> 0.2370744, 0.2624772, 3.0753158, 2.5379842, 0.6871498, 1.06…
$ Shape_Area <dbl> 0.0015239210, 0.0035311037, 0.3268678399, 0.0683785064, 0.0…
$ ADM2_EN    <chr> "Aba North", "Aba South", "Abadam", "Abaji", "Abak", "Abaka…
$ ADM2_PCODE <chr> "NG001001", "NG001002", "NG008001", "NG015001", "NG003001",…
$ ADM2_REF   <chr> "Aba North", "Aba South", "Abadam", "Abaji", "Abak", "Abaka…
$ ADM2ALT1EN <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
$ ADM2ALT2EN <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
$ ADM1_EN    <chr> "Abia", "Abia", "Borno", "Federal Capital Territory", "Akwa…
$ ADM1_PCODE <chr> "NG001", "NG001", "NG008", "NG015", "NG003", "NG011", "NG02…
$ ADM0_EN    <chr> "Nigeria", "Nigeria", "Nigeria", "Nigeria", "Nigeria", "Nig…
$ ADM0_PCODE <chr> "NG", "NG", "NG", "NG", "NG", "NG", "NG", "NG", "NG", "NG",…
$ date       <date> 2016-11-29, 2016-11-29, 2016-11-29, 2016-11-29, 2016-11-29…
$ validOn    <date> 2019-04-17, 2019-04-17, 2019-04-17, 2019-04-17, 2019-04-17…
$ validTo    <date> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ SD_EN      <chr> "Abia South", "Abia South", "Borno North", "Federal Capital…
$ SD_PCODE   <chr> "NG00103", "NG00103", "NG00802", "NG01501", "NG00302", "NG0…
$ geometry   <MULTIPOLYGON [°]> MULTIPOLYGON (((7.401109 5...., MULTIPOLYGON (…

From the attributed visible, we can see the HDE source (NGA) has a column called ‘ADM1_EN’ which can be used to filter for water points in Osun, Nigeria. However, this is not present in the geoBoundaries dataset. As a result, we will be using the Humanitarian Data Exchange source, where we get the files - nga_admbnda_adm2.

NGA sf data.frame consists of many redundent fields. The code chunk below uses select() of dplyr to retain column 3, 4, 8 and 9.

NGA <- NGA %>%
  select(c(3:4, 8:9))

We then use the filter() to filter out the polygon features of Osun.

NGA <- NGA %>% filter(ADM1_EN == "Osun")

Now, we use the st_crs() to check the coordinate system of the data. As we can see, it uses the WGS 84 coordinate system. The data is using a Geographic projected system, however, this is system is not appropriate since we need to use distance and area measures.

st_crs(NGA)
Coordinate Reference System:
  User input: WGS 84 
  wkt:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["latitude",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["longitude",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]

Hence, we use st_transform() and not st_set_crs() as st_set_crs() assigns the EPSG code to the dataframe, however, now we need to transform the dataframe from geographic to projected coordinate system. We will be using crs=26392 (found from the EPSGfor Nigeria).

NGA <- st_transform(NGA, crs = 26392)

Verify that the CRS of NGA dataframe has changed.

st_crs(NGA)
Coordinate Reference System:
  User input: EPSG:26392 
  wkt:
PROJCRS["Minna / Nigeria Mid Belt",
    BASEGEOGCRS["Minna",
        DATUM["Minna",
            ELLIPSOID["Clarke 1880 (RGS)",6378249.145,293.465,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4263]],
    CONVERSION["Nigeria Mid Belt",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",4,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",8.5,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.99975,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",670553.98,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Engineering survey, topographic mapping."],
        AREA["Nigeria between 6°30'E and 10°30'E, onshore and offshore shelf."],
        BBOX[3.57,6.5,13.53,10.51]],
    ID["EPSG",26392]]

3.2 Data Pre-processing

3.2.1 Dropping Invalid Dimensions

Since, we only have one dataframe, there are no invalid dimensions, and hence, this step is not required.

3.2.2 Invalid Geometries

The st_is_valid() function checks whether a geometry is valid and returns the indices. Whereas, the length gives you a count of the indices with invalid geometries.

length(which(st_is_valid(NGA) == FALSE))
[1] 0

None of the values are Invalid, so we are good to go!!

3.2.3 Checking for Duplicated Names

We need to check for duplicate name in the data main data fields. Using duplicated() of Base R, we can flag out LGA names that might be duplicated as shown in the code chunk below.

NGA$ADM2_EN[duplicated(NGA$ADM2_EN)==TRUE]
character(0)

There are no duplicated values, so we are good to go!

3.2.4 Initial Visualization

plot(st_geometry(NGA))

4.0 Data Wrangling : Aspatial Data

4.1 Importing Aspatial Data

Since the WPdx data is in CSV format, we will use read_csv() of readr package to import WPdx.csv. The output is called wp_nga and is a tibble dataframe

wp_nga <- read_csv("data/aspatial/WPdx.csv") %>%
  filter(`#clean_country_name` == "Nigeria" & `#clean_adm1` == "Osun")

4.2 Converting water point data into sf point features

Converting an aspatial data into an sf data.frame involves two steps.

First, we need to convert the wkt field into sfc field by using st_as_sfc() function. The function stores it in a tibble data format.

wp_nga$Geometry = st_as_sfc(wp_nga$`New Georeferenced Column`)
wp_nga
# A tibble: 5,557 × 71
   row_id `#source`      #lat_…¹ #lon_…² #repo…³ #stat…⁴ #wate…⁵ #wate…⁶ #wate…⁷
    <dbl> <chr>            <dbl>   <dbl> <chr>   <chr>   <chr>   <chr>   <chr>  
 1 429123 GRID3             8.02    5.06 08/29/… Unknown <NA>    <NA>    Tapsta…
 2  70566 Federal Minis…    7.32    4.79 05/11/… No      Protec… Well    Mechan…
 3  70578 Federal Minis…    7.76    4.56 05/11/… No      Boreho… Well    Mechan…
 4  66401 Federal Minis…    8.03    4.64 04/30/… No      Boreho… Well    Mechan…
 5 422190 GRID3             7.87    4.88 08/29/… Unknown <NA>    <NA>    Tapsta…
 6 422064 GRID3             7.7     4.89 08/29/… Unknown <NA>    <NA>    Tapsta…
 7  65607 Federal Minis…    7.89    4.71 05/12/… No      Boreho… Well    Mechan…
 8  68989 Federal Minis…    7.51    4.27 05/07/… No      Boreho… Well    <NA>   
 9  67708 Federal Minis…    7.48    4.35 04/29/… Yes     Boreho… Well    Mechan…
10  66419 Federal Minis…    7.63    4.50 05/08/… Yes     Boreho… Well    Hand P…
# … with 5,547 more rows, 62 more variables: `#water_tech_category` <chr>,
#   `#facility_type` <chr>, `#clean_country_name` <chr>, `#clean_adm1` <chr>,
#   `#clean_adm2` <chr>, `#clean_adm3` <chr>, `#clean_adm4` <chr>,
#   `#install_year` <dbl>, `#installer` <chr>, `#rehab_year` <lgl>,
#   `#rehabilitator` <lgl>, `#management_clean` <chr>, `#status_clean` <chr>,
#   `#pay` <chr>, `#fecal_coliform_presence` <chr>,
#   `#fecal_coliform_value` <dbl>, `#subjective_quality` <chr>, …

Next, we use the st_sf() to convert the tibble data.frame into an sf object. It is also important for us to include the referencing system of the data into the sf object.

wp_sf <- st_sf(wp_nga, crs=4326)
wp_sf
Simple feature collection with 5557 features and 70 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 4.032004 ymin: 7.060309 xmax: 5.06 ymax: 8.061898
Geodetic CRS:  WGS 84
# A tibble: 5,557 × 71
   row_id `#source`      #lat_…¹ #lon_…² #repo…³ #stat…⁴ #wate…⁵ #wate…⁶ #wate…⁷
 *  <dbl> <chr>            <dbl>   <dbl> <chr>   <chr>   <chr>   <chr>   <chr>  
 1 429123 GRID3             8.02    5.06 08/29/… Unknown <NA>    <NA>    Tapsta…
 2  70566 Federal Minis…    7.32    4.79 05/11/… No      Protec… Well    Mechan…
 3  70578 Federal Minis…    7.76    4.56 05/11/… No      Boreho… Well    Mechan…
 4  66401 Federal Minis…    8.03    4.64 04/30/… No      Boreho… Well    Mechan…
 5 422190 GRID3             7.87    4.88 08/29/… Unknown <NA>    <NA>    Tapsta…
 6 422064 GRID3             7.7     4.89 08/29/… Unknown <NA>    <NA>    Tapsta…
 7  65607 Federal Minis…    7.89    4.71 05/12/… No      Boreho… Well    Mechan…
 8  68989 Federal Minis…    7.51    4.27 05/07/… No      Boreho… Well    <NA>   
 9  67708 Federal Minis…    7.48    4.35 04/29/… Yes     Boreho… Well    Mechan…
10  66419 Federal Minis…    7.63    4.50 05/08/… Yes     Boreho… Well    Hand P…
# … with 5,547 more rows, 62 more variables: `#water_tech_category` <chr>,
#   `#facility_type` <chr>, `#clean_country_name` <chr>, `#clean_adm1` <chr>,
#   `#clean_adm2` <chr>, `#clean_adm3` <chr>, `#clean_adm4` <chr>,
#   `#install_year` <dbl>, `#installer` <chr>, `#rehab_year` <lgl>,
#   `#rehabilitator` <lgl>, `#management_clean` <chr>, `#status_clean` <chr>,
#   `#pay` <chr>, `#fecal_coliform_presence` <chr>,
#   `#fecal_coliform_value` <dbl>, `#subjective_quality` <chr>, …

Like step 3.2 Data Pre-processing, we transform the projection from wgs84 to the appropriate projected coordinate system of Nigeria.

wp_sf <- wp_sf %>%
  st_transform(crs = 26392)

4.3 Data Wrangling for Water Data Point

Exploratory Data Analysis (EDA) helps to gain initial understanding of the data. The freq() of funModeling package is used to reveal the distribution of water point status visually.

freq(data = wp_sf,
     input = '#status_clean')

                     #status_clean frequency percentage cumulative_perc
1                       Functional      2319      41.73           41.73
2                   Non-Functional      2008      36.13           77.86
3                             <NA>       748      13.46           91.32
4      Functional but needs repair       248       4.46           95.78
5 Non-Functional due to dry season       151       2.72           98.50
6        Functional but not in use        63       1.13           99.63
7                        Abandoned        15       0.27           99.90
8         Abandoned/Decommissioned         5       0.09          100.00

The diagram shows that there are nine classes present in the ‘status_clean’ field. Hence, now we will be performing data wrangling tasks to create 3 data object - Functional, Non-Functional and Unknown.

We use rename() function from the dplyr package to rename the column from #status_clean to status_clean for easier handling in subsequent steps. select() is used to include status_clean in the output sf data.frame. We use the mutate() and replace_na() functions to recode all the NA values in status_clean into unknown.

wp_sf_nga <- wp_sf %>% 
  rename(status_clean = '#status_clean') %>%
  select(status_clean) %>%
  mutate(status_clean = replace_na(
    status_clean, "unknown"))

4.3.1 Extracting Water Point Data

Now we are ready to extract the water point data according to their status.

The code chunk below is used to extract functional water point.

wp_functional <- wp_sf_nga %>%
  filter(status_clean %in%
           c("Functional",
             "Functional but not in use",
             "Functional but needs repair"))

The code chunk below is used to extract nonfunctional water point.

wp_nonfunctional <- wp_sf_nga %>%
  filter(status_clean %in%
           c("Abandoned/Decommissioned",
             "Abandoned",
             "Non-Functional due to dry season",
             "Non-Functional",
             "Non functional due to dry season"))

The code chunk below is used to extract water point with unknown status.

wp_unknown <- wp_sf_nga %>%
  filter(status_clean == "unknown")

Performing a quick EDA on the derived sfa.dataframes

freq(data = wp_functional,
     input = 'status_clean')

                 status_clean frequency percentage cumulative_perc
1                  Functional      2319      88.17           88.17
2 Functional but needs repair       248       9.43           97.60
3   Functional but not in use        63       2.40          100.00
freq(data = wp_nonfunctional,
     input = 'status_clean')

                      status_clean frequency percentage cumulative_perc
1                   Non-Functional      2008      92.15           92.15
2 Non-Functional due to dry season       151       6.93           99.08
3                        Abandoned        15       0.69           99.77
4         Abandoned/Decommissioned         5       0.23          100.00
freq(data = wp_unknown,
     input = 'status_clean')

  status_clean frequency percentage cumulative_perc
1      unknown       748        100             100

We can see from the map below, the proportion of functional and non-functional water is quite similar.

tmap_mode("view")
tm_shape(wp_functional) +
 tm_dots(col = "status_clean",
         pal = "blue",
         title = "Functional") +
tm_shape(wp_nonfunctional) +
 tm_dots(col = "status_clean",
         pal = "red",
         title = "Non-Functional") +
  tm_view(set.zoom.limits = c(8.5,15))

4.3.2 Performing Point-In Polygon Count

Next, we want to find out the number of total, functional, nonfunctional and unknown water points in Osun State. This is performed in the following code chunk. First, it identifies the functional water points in each LGA by using st_intersects() of sf package. Next, length() is used to calculate the number of functional water points that fall inside each LGA.

NGA_wp <- NGA %>% 
  mutate(`total_wp` = lengths(
    st_intersects(NGA, wp_sf_nga))) %>%
  mutate(`wp_functional` = lengths(
    st_intersects(NGA, wp_functional))) %>%
  mutate(`wp_nonfunctional` = lengths(
    st_intersects(NGA, wp_nonfunctional))) %>%
  mutate(`wp_unknown` = lengths(
    st_intersects(NGA, wp_unknown)))
NGA_wp
Simple feature collection with 30 features and 8 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 176503.2 ymin: 331434.7 xmax: 291043.8 ymax: 454520.1
Projected CRS: Minna / Nigeria Mid Belt
First 10 features:
          ADM2_EN ADM2_PCODE ADM1_EN ADM1_PCODE                       geometry
1        Aiyedade   NG030001    Osun      NG030 MULTIPOLYGON (((213526.6 34...
2        Aiyedire   NG030002    Osun      NG030 MULTIPOLYGON (((212542.6 40...
3  Atakumosa East   NG030003    Osun      NG030 MULTIPOLYGON (((265746.8 37...
4  Atakumosa West   NG030004    Osun      NG030 MULTIPOLYGON (((248871.4 40...
5      Boluwaduro   NG030005    Osun      NG030 MULTIPOLYGON (((266092.2 43...
6          Boripe   NG030006    Osun      NG030 MULTIPOLYGON (((255072.5 43...
7       Ede North   NG030007    Osun      NG030 MULTIPOLYGON (((236386.9 41...
8       Ede South   NG030008    Osun      NG030 MULTIPOLYGON (((236386.9 41...
9        Egbedore   NG030009    Osun      NG030 MULTIPOLYGON (((220756 4317...
10         Ejigbo   NG030010    Osun      NG030 MULTIPOLYGON (((214422.1 42...
   total_wp wp_functional wp_nonfunctional wp_unknown
1       389           157              154         78
2       175            89               57         29
3       223            98               92         33
4       246           111              103         32
5       129            63               51         15
6       177            79               85         13
7       216           141               50         25
8       146            72               39         35
9       142            63               44         35
10      434           274              126         34

We then visualise attributes by using statistcal graph. We use functions of ggplot2 package to reveal the distribution of total water points in Osun’s LGA using histogram.

ggplot(data = NGA_wp,
       aes(x = total_wp)) + 
  geom_histogram(bins=20,
                 color="black",
                 fill="light blue") +
  geom_vline(aes(xintercept=mean(
    total_wp, na.rm=T)),
             color="red", 
             linetype="dashed", 
             size=0.8) +
  ggtitle("Distribution of total water points") +
  xlab("No. of water points") +
  ylab("No. of\nLGAs") +
  theme(axis.title.y=element_text(angle = 0))

5.0 Combined Data Wrangling : Geospatial & Aspatial Data

This is an essential step, which is the process of getting data from its raw input into a format which can be used for anlaysis.

5.1 Converting sf data frames to sp’s Spatial* Class

We use the as_spatial() function to convert the three geospatial data from simple feature data frame to sp’s Spatial* class.

wp_functional_spatial = as_Spatial(wp_functional)
wp_nonfunctional_spatial = as_Spatial(wp_nonfunctional)
NGA_spatial <- as_Spatial(NGA)
NGA_spatial
class       : SpatialPolygonsDataFrame 
features    : 30 
extent      : 176503.2, 291043.8, 331434.7, 454520.1  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=4 +lon_0=8.5 +k=0.99975 +x_0=670553.98 +y_0=0 +a=6378249.145 +rf=293.465 +towgs84=-92,-93,122,0,0,0,0 +units=m +no_defs 
variables   : 4
names       :  ADM2_EN, ADM2_PCODE, ADM1_EN, ADM1_PCODE 
min values  : Aiyedade,   NG030001,    Osun,      NG030 
max values  :   Osogbo,   NG030030,    Osun,      NG030 
wp_functional_spatial
class       : SpatialPointsDataFrame 
features    : 2630 
extent      : 177285.9, 290751, 343128.1, 450859.7  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=4 +lon_0=8.5 +k=0.99975 +x_0=670553.98 +y_0=0 +a=6378249.145 +rf=293.465 +towgs84=-92,-93,122,0,0,0,0 +units=m +no_defs 
variables   : 1
names       :              status_clean 
min values  :                Functional 
max values  : Functional but not in use 
wp_nonfunctional_spatial
class       : SpatialPointsDataFrame 
features    : 2179 
extent      : 180539, 290616, 340054.1, 450780.1  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=4 +lon_0=8.5 +k=0.99975 +x_0=670553.98 +y_0=0 +a=6378249.145 +rf=293.465 +towgs84=-92,-93,122,0,0,0,0 +units=m +no_defs 
variables   : 1
names       :                     status_clean 
min values  :                        Abandoned 
max values  : Non-Functional due to dry season

5.2 Converting from Spatial* classes to sp format

In order to use the spatstat for our analysis, we need our data to be in the ppp object form. Hence, we first need to convert them into Spatial object first and then into ppp object.

# convert into respective sp (in our case, either polygons or points)
wp_functional_sp <- as(wp_functional_spatial, "SpatialPoints")
wp_nonfunctional_sp <- as(wp_nonfunctional_spatial, "SpatialPoints")
NGA_sp <-as(NGA_spatial, "SpatialPolygons")
wp_functional_sp
class       : SpatialPoints 
features    : 2630 
extent      : 177285.9, 290751, 343128.1, 450859.7  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=4 +lon_0=8.5 +k=0.99975 +x_0=670553.98 +y_0=0 +a=6378249.145 +rf=293.465 +towgs84=-92,-93,122,0,0,0,0 +units=m +no_defs 
wp_nonfunctional_sp
class       : SpatialPoints 
features    : 2179 
extent      : 180539, 290616, 340054.1, 450780.1  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=4 +lon_0=8.5 +k=0.99975 +x_0=670553.98 +y_0=0 +a=6378249.145 +rf=293.465 +towgs84=-92,-93,122,0,0,0,0 +units=m +no_defs 
NGA_sp
class       : SpatialPolygons 
features    : 30 
extent      : 176503.2, 291043.8, 331434.7, 454520.1  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=4 +lon_0=8.5 +k=0.99975 +x_0=670553.98 +y_0=0 +a=6378249.145 +rf=293.465 +towgs84=-92,-93,122,0,0,0,0 +units=m +no_defs

5.3 Converting from sp format to spatstat ppp format

We can’t convert SpatialPolygons to ppp format - nor is there any need to. Hence, we won’t be including our ‘base map’, NGA.

# from sp object, convert into ppp format
wp_functional_ppp <- as(wp_functional_sp, "ppp")
wp_nonfunctional_ppp <- as(wp_nonfunctional_sp, "ppp")

The below map shows the point paterns for both functional and non-functional water points.

par(mfrow=c(1,2))
plot(wp_nonfunctional_ppp)
plot(wp_functional_ppp)

5.3.1 Handling Duplicated Points + Jittering

any(duplicated(wp_functional_ppp)) 
[1] FALSE
any(duplicated(wp_nonfunctional_ppp)) 
[1] FALSE

Since there is no duplication, we dont have to apply the process of Jittering.

5.4 Creating Owin Object

We need to now confine the analysis with a geographical area - Osun State and we do this by creating a object called owin which represent the polygonal region. The below code covert the SpatialPolygon (NGA_sp) created into an owin object.

NGA_owin <- as(NGA_sp, "owin")
plot(NGA_owin)

5.5 Combining point events object and owin object

In this step, we extract the functional and non-functional water points that are located within Osun, Nigeria. This combines both the point and polygon feature into one ppp object class.

wp_functional_ppp = wp_functional_ppp[NGA_owin]
wp_nonfunctional_ppp = wp_nonfunctional_ppp[NGA_owin]
par(mfrow=c(1,2))
plot(wp_nonfunctional_ppp)
plot(wp_functional_ppp)

6.0 Exploratory Spatial Data Analysis

In this section, we use the Hands-on Excercise 04 to help

  • Derive kernel density maps of functional and non-functional water points. Using appropriate tmap functions,

  • Display the kernel density maps on openstreetmap of Osub State, Nigeria.

  • Describe the spatial patterns revealed by the kernel density maps. Highlight the advantage of kernel density map over point map.

6.1 Kernel Density Estimation

6.1.1 Computing Kernel Density Estimation

There are two types of bandwidth methods - Fixed (Automatic) and Adaptive bandwidth method. These methods employ different uniform bases in density calculation.

Computing using Automatic Bandwidth selection method

We are using the below code to compute the kernel density by using the bw.diggle() - an automatic selection method and the smoothing method - kernel.

kde_wpfunctional_bw <- density(wp_functional_ppp,
                              sigma=bw.diggle,
                              edge=TRUE,
                            kernel="gaussian") 
kde_wpnonfunctional_bw <- density(wp_nonfunctional_ppp,
                              sigma=bw.diggle,
                              edge=TRUE,
                            kernel="gaussian") 

par(mfrow=c(1,2))
plot(kde_wpfunctional_bw,
     main = "Functional Water Points",
     ribside=c("right"))
plot(kde_wpnonfunctional_bw,
     main = "Non-Functional Water Points",
     ribside=c("right"))

Computing using Adaptive Bandwidth selection method

kde_wpfunctional_adaptive <- adaptive.density(wp_functional_ppp, method="kernel")

kde_wpnonfunctional_adaptive <- adaptive.density(wp_nonfunctional_ppp, method="kernel")

par(mfrow=c(1,2))
plot(kde_wpfunctional_adaptive,
     main = "Functional Water Points",
     ribside=c("right"))
plot(kde_wpnonfunctional_adaptive,
     main = "Non-Functional Water Points",
     ribside=c("right"))

Comparing Automated and Adapting Bandwidth Methods (side-by-side)

par(mfrow=c(1,2))
plot(kde_wpfunctional_bw,
     main = "Functional Water Points - Automated",
     ribside=c("right"))
plot(kde_wpfunctional_adaptive,
     main = "Functional Water Points - Adaptive",
     ribside=c("right"))

par(mfrow=c(1,2))
plot(kde_wpnonfunctional_bw,
     main = "Non-Functional Water Points - Automated",
     ribside=c("right"))
plot(kde_wpnonfunctional_adaptive,
     main = "Non-Functional Water Points - Adaptive",
     ribside=c("right"))

6.1.2 Rescalling KDE Values

As we can the KDE values are small (ranging from 0 to 0.000035). This is because the default unit of measurement of svy21 is in meter. As a result, the density values computed is in “number of points per square meter”. So rescale() is used to covert the unit of measurement from meter to kilometer.

wp_functional_ppp_km <- rescale(wp_functional_ppp, 1000, "km")
wp_nonfunctional_ppp_km <- rescale(wp_nonfunctional_ppp, 1000, "km")

Now we re-plot the graphs

kde_wpfunctional_km <- density(wp_functional_ppp_km,
                              sigma=bw.diggle,
                              edge=TRUE,
                            kernel="gaussian") 
kde_wpnonfunctional_km <- density(wp_nonfunctional_ppp_km,
                              sigma=bw.diggle,
                              edge=TRUE,
                            kernel="gaussian") 

par(mfrow=c(1,2))
plot(kde_wpfunctional_bw,
     main = "Functional Water Points",
     ribside=c("right"))
plot(kde_wpnonfunctional_bw,
     main = "Non-Functional Water Points",
     ribside=c("right"))

kde_wpfunctional_adaptive_km <- adaptive.density(wp_functional_ppp_km, method="kernel")

kde_wpnonfunctional_adaptive_km <- adaptive.density(wp_nonfunctional_ppp_km, method="kernel")

par(mfrow=c(1,2))
plot(kde_wpfunctional_adaptive,
     main = "Functional Water Points",
     ribside=c("right"))
plot(kde_wpnonfunctional_adaptive,
     main = "Non-Functional Water Points",
     ribside=c("right"))

For this assignment, we will be using the Automated Bandwidth method because it defines its base in geographical space, where as the Adaptive method defines it in population (Reference). As we learned in Chapter 04, Automated Bandwidth is very sensitive to highly skew distribution of spatial point patterns over geographical units (e.g - urban versus rural). However, since we don’t have highly skewed data (as seen in the Distribution Graph of Water Points), we can use Fixed/Automated Bandwidth method.

6.2 Converting KDE output into grid object

gridded_wpfunctional <- as.SpatialGridDataFrame.im(kde_wpfunctional_km)
gridded_wpnonfunctional <- as.SpatialGridDataFrame.im(kde_wpnonfunctional_km)

spplot(gridded_wpfunctional)

spplot(gridded_wpnonfunctional)

6.2.1 Converting Gridded Output into Raster

kde_wpfunctional_raster <- raster(gridded_wpfunctional)
kde_wpfunctional_raster
class      : RasterLayer 
dimensions : 128, 128, 16384  (nrow, ncol, ncell)
resolution : 0.8948485, 0.9616045  (x, y)
extent     : 176.5032, 291.0438, 331.4347, 454.5201  (xmin, xmax, ymin, ymax)
crs        : NA 
source     : memory
names      : v 
values     : -5.092436e-15, 25.49435  (min, max)
kde_wpnonfunctional_raster <- raster(gridded_wpnonfunctional)
kde_wpnonfunctional_raster
class      : RasterLayer 
dimensions : 128, 128, 16384  (nrow, ncol, ncell)
resolution : 0.8948485, 0.9616045  (x, y)
extent     : 176.5032, 291.0438, 331.4347, 454.5201  (xmin, xmax, ymin, ymax)
crs        : NA 
source     : memory
names      : v 
values     : -3.925434e-15, 20.49412  (min, max)

6.2.2 Assigning Projection Systems

projection(kde_wpfunctional_raster) <- CRS("+init=EPSG:26392 +datum:WGS84 +units=km")
kde_wpfunctional_raster
class      : RasterLayer 
dimensions : 128, 128, 16384  (nrow, ncol, ncell)
resolution : 0.8948485, 0.9616045  (x, y)
extent     : 176.5032, 291.0438, 331.4347, 454.5201  (xmin, xmax, ymin, ymax)
crs        : +proj=tmerc +lat_0=4 +lon_0=8.5 +k=0.99975 +x_0=670553.98 +y_0=0 +a=6378249.145 +rf=293.465 +units=km +no_defs 
source     : memory
names      : v 
values     : -5.092436e-15, 25.49435  (min, max)
projection(kde_wpnonfunctional_raster) <- CRS("+init=EPSG:26392 +datum:WGS84 +units=km")
kde_wpnonfunctional_raster
class      : RasterLayer 
dimensions : 128, 128, 16384  (nrow, ncol, ncell)
resolution : 0.8948485, 0.9616045  (x, y)
extent     : 176.5032, 291.0438, 331.4347, 454.5201  (xmin, xmax, ymin, ymax)
crs        : +proj=tmerc +lat_0=4 +lon_0=8.5 +k=0.99975 +x_0=670553.98 +y_0=0 +a=6378249.145 +rf=293.465 +units=km +no_defs 
source     : memory
names      : v 
values     : -3.925434e-15, 20.49412  (min, max)

6.3 Kernel Density Maps on OpenStreetMap

Now, as the assignment requirements has specified, we should plot our kernel density maps on OpenStreetMap. Since we’ll be plotting a lot of kernel density maps, let’s create a function:

density_map <- function(raster_object, map_title) {
  tmap_mode("view")
  tm_basemap("OpenStreetMap") +
tm_shape(raster_object) +
  tm_raster("v", alpha=0.9) + 
  tm_layout(legend.position = c("right", "bottom"), 
            legend.height = 0.5, 
            legend.width = 0.4,
            main.title = map_title,
            main.title.position = 'center',
            main.title.size = 1,
            frame = TRUE) + 
  tm_view(set.zoom.limits = c(8, 13))
  } 
kde_wpfunctional_density_map <- density_map(kde_wpfunctional_raster, map_title = "Functional Water Points in Osun State")
kde_wpnonfunctional_density_map <- density_map(kde_wpnonfunctional_raster, map_title = "Non-Functional Water Points in Osun State")
kde_wpfunctional_density_map
kde_wpnonfunctional_density_map
tmap_mode('plot')
tm_basemap("OpenStreetMap") +
tm_shape(kde_wpfunctional_raster) +
  tm_raster("v")

tmap_mode('plot')
tm_basemap("OpenStreetMap") +
tm_shape(kde_wpnonfunctional_raster) +
  tm_raster("v")

6.4 Kernel Density Maps Analysis

As we can see in the map in [5.5 Combining point events object and owin object], both the plots are similar with the Functional Water Point being comparatively a bit more denser (more point) than the Non-Functional Water Point. From the maps above, we can see that both the Functional and Non-Functional waterpoints are spread out, however, they are more concentrated in the center and the upper part of Osun. We don’t see that many waterpoints in lower part of Osun.

The Functional Water Points are slightly more spread out compared to the Non Functional Water Points, however, what is interesting to note is that the points in both the maps kind of coincide with each other. That is the points in the Functional Water Point seem to be close to that of the Non-Functional Water Point.

6.5 Advantage of Kernel Density Map over Point Map

To understand the advantage of Kernel Density Map over Point Map, we first need to plot the two and compare the differences.

tmap_mode("plot")
tm_shape(NGA_wp) +
  tm_borders(alpha = 0.5) +
  tmap_options(check.and.fix = TRUE) +
tm_shape(wp_nonfunctional) +
  tm_dots(col="red", size=0.05) +
  tm_layout(main.title = "Non-Functional Water Points",
          main.title.position = "center",
          main.title.size = 1.2,
          frame = TRUE)

kde_wpnonfunctional_density_map

With the Kernel Density Map, denser areas with a heavier distribution of Non-Functional Water Points are easily spotted. This is because the kernel density z-estimate helps to smooth out the points in a given area. Compared to the point map which just shows the points. Further, the gradient colour available (ranging from yellow to green) helps in understanding the density/concentration of water pumps in the area. It clearly shows the viewer which are the areas with more non-functional water pumps, however, with the point map, the users have to gauge/estimate which are the densers with more non-functional water points.

Hence to conclude, the Kernal Density provides a quantitative value representing the concentration of points, where as this can only be observed/gauged in Point Map.

  • With kernel density maps, it takes into account the inverse-distance-weighted counts of points, to represent the concentration of points at a particular location. This cannot be achieved through observation using point maps.

  • The inverse-distance-weighted counts is important because in the real-world, childcare centres that are further away from a particular location does not mean that they cannot potentially serve the population. These points should still be taken into account, just that points further away should just be given less weight, as people will have to travel further to access the childcare service. This is exactly what is accounted for with kernel function.

6.6?? Nearest Neighbour Analysis

The 95% confident interval will be used.

The test hypotheses for Functional Water Point is :

H0 : The distribution of Functional Water Point in Osun State is randomly distributed.

H1 : The distribution of Functional Water Point in Osun State is not randomly distributed.

clarkevans.test(wp_functional_ppp,
                correction="none",
                clipregion="NGA_owin",
                alternative=c("clustered"),
                nsim=99)

    Clark-Evans test
    No edge correction
    Monte Carlo test based on 99 simulations of CSR with fixed n

data:  wp_functional_ppp
R = 0.44265, p-value = 0.01
alternative hypothesis: clustered (R < 1)

Conclusion :

The test hypotheses for Non-Functional Water Point is :

H0 : The distribution of Non-Functional Water Point in Osun State is randomly distributed.

H1 : The distribution of Non-Functional Water Point in Osun State is not randomly distributed.

clarkevans.test(wp_nonfunctional_ppp,
                correction="none",
                clipregion="nga_owin",
                alternative=c("clustered"),
                nsim=99)

    Clark-Evans test
    No edge correction
    Monte Carlo test based on 99 simulations of CSR with fixed n

data:  wp_nonfunctional_ppp
R = 0.43223, p-value = 0.01
alternative hypothesis: clustered (R < 1)

Conclusion :

6.7 Colocation of Functional and Non-Functional Water Points

wp_sf_withoutUnknown <- wp_sf_nga %>%  filter(!status_clean=='unknown')
# This is required for Take Home Excercise 3
nb <- include_self(
  st_knn(st_geometry(wp_sf_withoutUnknown),6))

wt <- st_kernel_weights(nb,
                        wp_sf_withoutUnknown,
                        "gaussian",
                        adaptive = TRUE)

A <- wp_functional$status_clean

B <- wp_nonfunctional$status_clean

LCLQ <- local_colocation(A, B, nb, wt, 49)

LCLQ_stores <- cbind(wp_sf_withoutUnknown, LCLQ)
tmap_mode("view")
tm_shape(NGA) + 
  tm_polygons() + 
  tm_shape(LCLQ_stores)+
    tm_dots(col = "Non.Functional",
            size = 0.05,
            border.col = "grey",
            border.lwd = 0.5) + 
  tm_view(set.zoom.limits = c(8,11))
tmap_mode("view")
tm_shape(NGA) + 
  tm_polygons() + 
  tm_shape(LCLQ_stores)+
    tm_dots(col = "Abandoned",
            size = 0.05,
            border.col = "grey",
            border.lwd = 0.5) + 
  tm_view(set.zoom.limits = c(9,13))
tmap_mode("view")
tm_shape(NGA) + 
  tm_polygons() + 
  tm_shape(LCLQ_stores)+
    tm_dots(col = "Non.Functional.due.to.dry.season",
            size = 0.05,
            border.col = "grey",
            border.lwd = 0.5) + 
  tm_view(set.zoom.limits = c(9,11))

The above maps show the colocation of Functional Water Points and Non-Functional Water Points. ELABORATEEE MOREE

7.0 Second-order Spatial Point Patterns Analysis

For Functional Water Point in Osun Sate

  • H0: The distribution of the Functional Water Points are randomly distributed

  • H1: The distribution of the Functional Water Points are not randomly distributed

  • Confidence level : 99%

  • Significance level : 0.01 In light that 0.05 is the most common level of significance, I’ve decided to make it slightly stricter and use 0.01 - I believe a 1% risk of an incorrect hypothesis is good enough, especially considering the trade-offs between sensitivity and false positives for this hypothesis testing

For Non-Functional Water Point in Osun Sate

  • H0: The distribution of the Non-Functional Water Points are randomly distributed

  • H1: The distribution of the Non-Functional Water Points are not randomly distributed

  • Confidence level : 99%

  • Significance level : 0.01 In light that 0.05 is the most common level of significance, I’ve decided to make it slightly stricter and use 0.01 - I believe a 1% risk of an incorrect hypothesis is good enough, especially considering the trade-offs between sensitivity and false positives for this hypothesis testing

7.1 Analysing Spatial Point Process Using G-Function

7.1.1 Functional Water Point

Computing G-function estimation

G_wp_functional = Gest(wp_functional_ppp, correction = "border")
plot(G_wp_functional, xlim=c(0,500))

Performing Complete Spatial Randomness Test

G_wp_functional.csr <- envelope(wp_functional_ppp, Gest, nsim = 999)
Generating 999 simulations of CSR  ...
1, 2, 3, ......10.........20.........30.........40.........50.........60........
.70.........80.........90.........100.........110.........120.........130......
...140.........150.........160.........170.........180.........190.........200....
.....210.........220.........230.........240.........250.........260.........270..
.......280.........290.........300.........310.........320.........330.........340
.........350.........360.........370.........380.........390.........400........
.410.........420.........430.........440.........450.........460.........470......
...480.........490.........500.........510.........520.........530.........540....
.....550.........560.........570.........580.........590.........600.........610..
.......620.........630.........640.........650.........660.........670.........680
.........690.........700.........710.........720.........730.........740........
.750.........760.........770.........780.........790.........800.........810......
...820.........830.........840.........850.........860.........870.........880....
.....890.........900.........910.........920.........930.........940.........950..
.......960.........970.........980.........990........ 999.

Done.
plot(G_wp_functional.csr)

Conclusion: The observed G(r) is far above the G(theo) as well as the envelope - indicating that Functional Water Points are clustered. Hence, we reject the null hypothesis that Functional Water Points are randomly distributed at 99% confident interval.

7.1.2 Non-Functional Water Point

Computing G-function estimation

G_wp_nonfunctional = Gest(wp_nonfunctional_ppp, correction = "border")
plot(G_wp_nonfunctional, xlim=c(0,500))

Performing Complete Spatial Randomness Test

G_wp_nonfunctional.csr <- envelope(wp_functional_ppp, Gest, nsim = 999)
Generating 999 simulations of CSR  ...
1, 2, 3, ......10.........20.........30.........40.........50.........60........
.70.........80.........90.........100.........110.........120.........130......
...140.........150.........160.........170.........180.........190.........200....
.....210.........220.........230.........240.........250.........260.........270..
.......280.........290.........300.........310.........320.........330.........340
.........350.........360.........370.........380.........390.........400........
.410.........420.........430.........440.........450.........460.........470......
...480.........490.........500.........510.........520.........530.........540....
.....550.........560.........570.........580.........590.........600.........610..
.......620.........630.........640.........650.........660.........670.........680
.........690.........700.........710.........720.........730.........740........
.750.........760.........770.........780.........790.........800.........810......
...820.........830.........840.........850.........860.........870.........880....
.....890.........900.........910.........920.........930.........940.........950..
.......960.........970.........980.........990........ 999.

Done.
plot(G_wp_nonfunctional.csr)

Conclusion: The observed G(r) is far above the G(theo) as well as the envelope - indicating that Non Functional Water Points are clustered. Hence, we reject the null hypothesis that Non Functional Water Points are randomly distributed at 99% confident interval.

7.2 Analysing Spatial Point Process Using L-Function

7.2.1 Functional Water Point

#|eval: false
#L_wp = Lest(wp_functional_ppp, correction = "Ripley")
#plot(L_wp, . -r ~ r, 
     #ylab= "L(d)-r", xlab = "d(m)")
#|eval: false
#L_wp.csr <- envelope(wp_functional_ppp, Lest, nsim = 39, rank = 1, glocal=TRUE)
#|eval: false
#plot(L_wp.csr, . - r ~ r, xlab="d", ylab="L(d)-r")

Non functional L test

#L_nonwp = Lest(wp_nonfunctional_ppp, correction = "Ripley")
#plot(L_ck_wp, . -r ~ r, 
     #ylab= "L(d)-r", xlab = "d(m)")
#|eval: false
#L_nonwp.csr <- envelope(wp_nonfunctional_ppp, Lest, nsim = 39, rank = 1, glocal=TRUE)
#|eval: false
#plot(L_nonwp.csr, . - r ~ r, xlab="d", ylab="L(d)-r")

8.0 Spatial Correlation Analysis

8.1 Data Pre-Processing

Convert sf data frames to sp’s Spatial class

#wp_spatial <- as_Spatial(wp_sf_nga)

Convert spatial class into generic sp class

#wp_sp <- as(wp_spatial, "SpatialPoints")

Converting generic sp format into spatstat’s ppp format

#wp_ppp <- as(wp_sp, "ppp")
#wp_ppp
#plot(wp_ppp)
#wp_ppp_km <- rescale(wp_ppp, 1000, "km")

9.0 References